home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / intro / dem01.dem
Text File  |  1999-09-16  |  5KB  |  191 lines
  1. mode(7)
  2. //               SCILAB OBJECTS
  3. //               1. SCALARS
  4. a=1               //constant
  5. 1==1              //boolean
  6. 'string'          //character string
  7. z=poly(0,'z')     // polynomial with variable 'z' and with one root at zero
  8. p=1+3*z+4.5*z^2   //polynomial 
  9. r=z/p             //rational
  10.  
  11. //                2. MATRICES
  12. a=[a+1 2 3
  13.      0 0 atan(1)
  14.      5 9 -1]      //constant matrix
  15.  
  16. b=[%t,%f]         //boolean matrix
  17.  
  18. mc=['this','is';
  19.     'a' ,'matrix']   //matrix of strings
  20.  
  21. mp=[p,1-z;
  22.     1,z*p]        //polynomial matrix
  23.  
  24. mp=[p 1-z]
  25. mp=[mp;1 1+z*p]   //matrix polynomial
  26.  
  27. f=mp/poly([1+%i 1-%i 1],'z')   //rational matrix
  28.  
  29. //                 3. LISTS
  30. l=list(a,-(1:5), mp,['this','is';'a','list'])   //list
  31. b=[1 0;0 1;0 0];c=[1 -1 0];d=0*c*b;x0=[0;0;0];
  32. sl=syslin('c',a,b,c,d,x0)    //Linear system in state-space representation.
  33. slt=ss2tf(sl)                // Transfer matrix
  34.  
  35. //                  OPERATIONS
  36. v=1:5;v*v'                   //constant matrix
  37. mp'*mp+eye                   //polynomial matrix
  38. mp1=mp(1,1)+4.5*%i           //complex
  39. fi=c*(z*eye-a)^(-1)*b;       //transfer function evaluation
  40. f(:,1)*fi                    //rationals
  41. m=[mp -mp; mp' mp+eye]       //usual Matlab syntax for polynomials
  42. [fi, fi(:,1)]                // ... or rationals
  43. f=syslin('c',f);             
  44. num=f(2);den=f(3);           //operation on transfer matrix
  45.  
  46. //                  SOME NUMERICAL PRIMITIVES
  47. inv(a)                       //Inverse
  48. inv(mp)                      //Inverse
  49. inv(sl*sl')                  //Product of two linear systems and inverse
  50. w=ss2tf(ans)                 //Transfer function representation
  51. inv(ss2tf(sl)*ss2tf(sl'))    //Product of two transfer functions and inverse
  52. clean(w-ans)                 
  53. n=contr(a,b)                 //Controllability
  54. k=ppol(a,b,[-1-%i -1+%i -1])        //Pole placement
  55. poly(a-b*k,'z')-poly([-1-%i -1+%i -1],'z')    //Check...
  56.  
  57. s=sin(0:0.1:5*%pi);
  58. ss=fft(s(1:128),-1);        //FFT
  59. xbasc();
  60. plot2d3("enn",1,abs(ss)');      //simple plot
  61.  
  62. x=lyap(a,diag([1 2 3]),'cont')   //Lyapunov equation
  63.  
  64.  
  65. //              ON LINE DEFINITION OF MACRO
  66. deff('[x]=fact(n)','if n=0 then x=1,else x=n*fact(n-1),end')
  67. 10+fact(5)
  68. //                    OPTIMIZATION
  69. deff('[f,g,ind]=rosenbro(x,ind)', 'a=x(2)-x(1)^2 , b=1-x(2) ,...
  70. f=100.*a^2 + b^2 , g(1)=-400.*x(1)*a , g(2)=200.*a -2.*b ');
  71. comp(rosenbro);[f,x,g]=optim(rosenbro,[2;2],'qn')
  72.  
  73. //                   SIMULATION
  74. a=rand(3,3)
  75. e=exp(a)
  76. deff('[ydot]=f(t,y)','ydot=a*y');comp(f)
  77. e(:,1)-ode([1;0;0],0,1,f)
  78.  
  79. //                  SYSTEM DEFINITION
  80. s=poly(0,'s')
  81. h=[1/s,1/(s+1);1/s/(s+1),1/(s+2)/(s+2)]
  82. w=tf2ss(h);
  83. ss2tf(w)
  84. h1=clean(ans)
  85.  
  86. //             EXAMPLE: SECOND ORDER SYSTEM ANALYSIS
  87. sl=syslin('c',1/(s*s+0.2*s+1))
  88. instants=0:0.05:20;
  89. //             step response:
  90. y=csim('step',instants,sl);
  91. xbasc();plot2d(instants',y')
  92. //             Delayed step response
  93. deff('[in]=u(t)','if t<3 then in=0;else in=1;end');
  94. comp(u);
  95. y1=csim(u,instants,sl);plot2d(instants',y1');
  96.  
  97. //             Impulse response;
  98. yi=csim('imp',instants,sl);xbasc();plot2d(instants',yi');
  99. yi1=csim('step',instants,s*sl);plot2d(instants',yi1');
  100.  
  101. //              Discretization
  102. dt=0.05;
  103. sld=dscr(tf2ss(sl),0.05);
  104.  
  105. //               Step response
  106. u=ones(instants);
  107. yyy=flts(u,sld);
  108. xbasc();plot(instants,yyy)
  109.  
  110. //              Impulse response
  111. u=0*ones(instants);u(1)=1/dt;
  112. yy=flts(u,sld);
  113. xbasc();plot(instants,yy)
  114.  
  115. //            system interconnexion
  116. w1=[w,w];
  117. clean(ss2tf(w1))
  118. w2=[w;w];
  119. clean(ss2tf(w2))
  120.  
  121. //               change of variable
  122. z=poly(0,'z');
  123. horner(h,(1-z)/(1+z))  //bilinear transform
  124.  
  125. //                 PRIMITIVES
  126. H=[1.    1.    1.    0.;
  127.    2.  - 1.    0.    1;
  128.    1.    0.    1.    1.;
  129.    0.    1.    2.  - 1];
  130.  
  131. ww=spec(H)
  132.  
  133. //             STABLE SUBSPACES
  134. [X,d]=schur(H,'cont');
  135. X'*H*X
  136.  
  137. [X,d]=schur(H,'disc');
  138. X'*H*X
  139.  
  140. //Selection of user-defined eigenvalues (# 3 and 4 here);
  141. deff('[flg]=sel(x)','flg=0,ev=x(2)/x(3),if abs(ev-ww(3))<0.0001|abs(ev-ww(4))<0.00001 then flg=1,end')
  142. [X,d]=schur(H,sel)
  143.  
  144. X'*H*X
  145.  
  146. //               With matrix pencil
  147. [X,d]=gschur(H,eye(H),sel)
  148. X'*H*X
  149.  
  150. //            block diagonalization
  151. [ab,x,bs]=bdiag(H);
  152.  
  153. inv(x)*H*x
  154.  
  155. //                     Matrix pencils
  156. E=rand(3,2)*rand(2,3);
  157. A=rand(3,2)*rand(2,3);
  158. s=poly(0,'s');
  159.  
  160. w=det(s*E-A)   //determinant
  161. [al,be]=gspec(A,E);
  162. al./(be+%eps*ones(be))
  163. roots(w)
  164. [Ns,d]=coffg(s*E-A);    //inverse of polynomial matrix;
  165. clean(Ns/d*(s*E-A))
  166. [Q,M,i1]=pencan(E,A);   // Canonical form;
  167. M*E*Q
  168. M*A*Q
  169.  
  170. //           PAUSE-RESUME
  171. write(%io(2),'pause command...');
  172. write(%io(2),'TO CONTINUE...');
  173. write(%io(2),'ENTER ''resume (or return) or click on resume!!''');
  174. //pause;
  175.  
  176. //           CALLING EXTERNAL ROUTINE
  177. foo=['      subroutine foo(a,b,c)';
  178.      '      c=a+b';
  179.      '      end'  ];
  180. unix_s('rm -f foo.f')
  181. write('foo.f',foo);
  182. unix_s('make foo.o')     //Compiling...(needs fortran compiler)
  183. //..................WARNING.......................
  184. //NEXT COMMAND LINE WILL DO THE LINK OF THE ROUTINE foo WITH SCILAB
  185. //THIS COMMAND MAY FAIL FOR SystemV COMPILED VERSIONS OF SCILAB 
  186. //BECAUSE LINK NEEDS THE LIBRARIES (SEE THE HELP OF "LINK" COMMAND)
  187. link('foo.o','foo')    //Linking to Scilab
  188. //5+7 by fortran
  189. fort('foo',5,1,'r',7,2,'r','out',[1,1],3,'r')
  190.  
  191.